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Abstract 

This paper consists of two parts. The first part proposes a new methodological framework 
within which the heat conductivity in ID lattices can be studied. The total process of heat 
QQ i conductivity is decomposed into two contributions where the first one is the equilibrium process 

at equal temperatures T of both lattice ends and the second - non-equilibrium process with the 
Q ■ temperature AT of one end and zero temperature of the other. This approach allows to isolate and 

u. 

^ . analyze the heat transfer in explicit form. The heat conductivity in the limit AT — t- is reduced to 

the heat conductivity of harmonic lattice with stochastic rigidities determined by the equilibrium 
process at temperature T. A threshold temperature Tthr is found which separates two regimes: 
t> . small perturbations exponentially decay at T < Tthr and tend to constant value at T > Tthr- 

m 

, The threshold temperature scales Tthr(-^) ~ A^^^ with the lattice size N and Tthr — in the 

' thermodynamic limit. Some unusual properties of heat conductivity can be exhibited on nanoscales 

' at low temperatures. The thermodynamics of the /3-FPU lattice can be adequately approximated 

by the harmonic lattice with temperature renormalized coefficients of rigidity. The second part 
testifies in the favor of the soliton and breather contribution to the heat conductivity in contrast 
to conclusions made in [N. Li, B. Li, S. Flach, PRL 105 (2010) 054102]. In the long-wavelength 
continuum limit the discrete /3-FPU lattice is reduced to the modified Korteweg - de Vries equation. 
This equation has soliton and breather solutions. Numerical simulations demonstrate their high 
stability. New method for the visualization of moving solitons and breathers is suggested. An 
accurate expression for the dependence of the sound velocity on temperature is also obtained. Our 
results support the conjecture on the solitons and breathers contribution to the heat conductivity. 
The fraction of total heat flux transferred by solitons and breathers merits additional analysis. 
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I. INTRODUCTION 



The problem of heat conductivity in low dimensional systems attracts much attention 
in last decades (see review and is motivated by the discovery of quasi-one-dimensional 
(nanotubes, nanowires, etc.) and two-dimensional (graphen, graphan, etc.) systems. 

The modern theory of heat conductivity was initiated by the celebrated preprint of 
E. Fermi, J. Pasta, and S. Ulam though the primary aim was "o/ establishing, experi- 
mentally, the rate of approaching to the equipartition of energy among the various degrees 
of freedom" . Subsequent investigations demonstrated wide area of consequences in many 
physical and mathematical phenomena (see reviews in special issues of journals CHAOS |3| 
and Lecture Notes in Physics j4| devoted to the 50th anniversary of the FPU preprint). 

The dynamical properties of nonlinear systems in microcanonical ensemble (total energy 
E = const) were thoroughly analyzed. It allows to investigate the dynamics and to get exact 
results (soliton js-?! and breather j8-12| solutions), to analyze regular and stochastic regimes 



and to find the corresponding thresholds. The 



in the field of "experimental mathematics" 
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'PU preprint also initiated the investigations 



About ten decades ago P. Debye argued that the nonlinearity can be responsible for the 
finite value of heat conductivity in insulating materials [l^ . But modern analysis shows that 
it is not always the case. There are many examples where the coefficient of heat conductivity 
X diverges with the increasing system size L as x oc L° where a > 0, and x — oo in the 
thermodynamic limit {L — )■ oo). Most of momentum conserving one-dimensional nonlinear 
lattices with various types of nearest-neighbor interactions have this unusual property (see, 
e:g., QQ, 



16| ). Moreover, some other systems, - two- 17H19| and three-dimensional lattices 
|20|, polyethylene chain (21], carbon nanotubes 22|-|26| have analogous property - diverging 
heat conductivity with the increasing size of the system. 

There were some conjectures explaining the anomalous heat conductivity. Generally 
speaking, whenever the equilibrium dynamics of a lattice can be decomposed into that of 
independent "modes" or quasi-particles, the system is expected to behave as an ideal thermal 
conductor j^^]. Thereby, the existence of stable nonlinear excitations is expected to yield 
ballistic rather than diffusive transport. At low temperatures normal modes are phonons. 
At higher temperatures noninteracting "gas" of solitons or/and breathers starts to play more 
significant role, and M. Toda was the first who suggested the possibility of heat transport 
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by solitons [28 1. 

Though analytical expressions for solitons can be derived only for few continuum models 
described by partial differential equations, Friesecke and Pego in a series of recent papers ^ 



32| made a detailed study of the existence and stability of solitary wave solutions on discrete 



lattices with the Hamiltonian H = + u{yi), where tji = Xi — pi = Xj. It has been 

proven that the systems with this Hamiltonian and with the following generic properties of 
nearest-neighbor interactions: u'{0) = 0; u"{0) > 0; u"'{0) ^ has a family of solitary wave 
solutions which in the small amplitude, long-wavelength limit have a profile close to that of 
the KdV soliton. It was also shown [33] that these solutions are asymptotically stable. Thus 
most acceptable point of view on the origin of anomalous heat conductivity in nonlinear 
lattices is as follows: phonons are responsible for heat conductivity at low temperatures. 



and at high temperatures - solitons 
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35|. 



Rather confusing experimental and numerical results are demonstrated in literature about 
the dependence of heat conductivity on different parameters, - model under consideration, 
types of boundary conditions, used thermostat and temperature. For instance, temperature 
dependence of heat conductivity in carbon nanotubes decreases as x ~ 1/T at T > 10 K 36| ; 
experimentally is found [23] that x also decreases with the growth of temperature. Different 
temperature dependencies x vs. T were found in ID ^nonlinear lattices. For /3-FPU lattice: 
X ~ iV°T^^ at T < 0.1 and x ~ N°'T^^^ at T > 50 [37] what is usually observed in insulating 
crystals. For the interparticle harmonic potentials and on-site potentials (e.g. Klein-Gordon 
chains) x ~ T^^'^^, i.e. heat conductivity decreases with the growth of temperature 38). 
But there exists firm theoretical background [39] that the exponent a in the dependence 
X oc A^" is the universal constant a 1/3 in momentum-conserving systems. 

The calculation of heat conductivity at small temperature gradients is an additional 
problem. Usually these calculations are very time consuming because of great fluctuations 
of heat current and statistical averaging over large number of MD trajectories is necessary. 

The paper organized as follows. In section [III the heat conductivity is considered when 
the temperature gradient VT is small. The explicit contribution to the heat conductivity 
is extracted by the decomposing of the total process into two parts. The first one is the 
equilibrium process at temperatures T of both lattice ends, and the second - non-equilibrium, 
when one lattice end has temperature AT and the other - zero temperature. Namely the 
latter process is responsible for the heat conductivity. This method allows to find the 
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threshold temperature Tthr- And though Tthr scales with the lattice length N as Tthr ~ A^~^, 
unusual dynamics can be revealed on nanoscale when both T and are small. The low 
temperature thermodynamics can be adequately described in terms of harmonic lattice with 
temperature renormalized rigidity coefficients. 

Section IIIII is independent of the previous one and is aimed at elucidating the role of 
solitons and breathers in the heat conductivity. Solitons and breathers are found as the 
solutions of the modified Korteweg - de Vries equation. This equation is obtained in the 
continuum limit from the discrete /3-FPU lattice. Soliton and breather solution are checked 
in numerical simulation and demonstrate very high stability. 

Necessary details of the derivation of accurate soliton and breather solutions in the con- 
tinuum limit are given in Appendix. 



II. HEAT CONDUCTIVITY IN THE /3-FPU LATTICE 

We consider the one- dimensional /3-FPU lattice of oscillators with the interaction of 
nearest neighbors 

U = '^'^{xi- Xi-iY + ^{xi- Xi-i)'^ (1) 

i 

(usually the dimensionless potential will be used below, that is a = (3 = m = 1). Nonequi- 
librium conditions are necessary for the heat transport simulation. The most abundant 
method is the placement of the lattice into the heat bath with different temperatures of left 
T+ and right T_ ends (T+ > T_). Different types of heat reservoirs are thoroughly analyzed 
in The usage of the Langevin forces with the noise terms and friction forces acting on 
the left F+ = ^+ — 7x1 and right F_ = ^-7Xiv oscillators is the common practice (7 = 1 is 
also put for brevity). {^±} are independent Wiener processes with zero mean and the cor- 
relator (^±(ti) '^±(^2)) = 2T± 6(ti — 12). AT = (T+ — T_) is the temperature difference. The 
generalized Langevin dynamics with a memory kernel and colored noises is also suggested 



40| to correctly account for the effect of the heat baths. 

The following set of stochastic differential equations (SDEs) 

dU 

x» = + SaF+ + 6iNF_ (2) 

is usually solved to find the heat flux J. And the local heat flux (power transmitted from 
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FIG. 1: Schematic representation of the total process x(t) as sum of equihbrium x^{t) and non- 
equiUbrium x^(t) processes. 



zth to {i + l)th oscillator) is 41 1 



(3) 



where Fi^i^i is a shorthand notation for the force exerted by the ith on the {i + l)th 
oscillator. The total heat flux J can be found as the mean J = {N — 1)^^ Xlf^"^ Ji-^i+i- 



A. Equilibrium and non-equilibrium contributions to the heat conductivity 

If T_ 7^ then the process of heat conductivity can be formally decomposed into two 
contributions: the first one - equilibrium process with equal temperatures T_ of both lattice 
ends; and the second - nonequilibrium process with temperature AT of the left lattice end 
and zero temperature of the right end (see Fig. [1]) (by 'process' we hereafter assume for 
brevity the solution x(t) = {xi(t), X2(t), . . . , x^it)} ; w{t) = {vi(t) , V2(t) , . . . , v^it)} of the 
corresponding SDKs). 

Namely the second process is responsible for the heat transport taking place on(?) the 
background of equilibrium process. Once this approach is utilized then the noise terms in ([21), 
owing to their independence, are {^+} = {^°} + {^^} and {^_} = {^°} for the left and right 
lattice ends, correspondingly Superscripts '0' and '1' refer to equilibrium and nonequilibrium 
processes. The total process x(t) can be represented as the sum 

x(t)=x°(t)+xi(t), (4) 

where x'^(t) is the equilibrium (Gibbs's) process at temperature T_, and x^(t) - nonequilib- 
rium, responsible for the energy transport, process. The corresponding stochastic dynamics 



IS 



X 



+5a(e°-i?) + 5.7v(e°-4), (5) 



dU dU' 



Q^o -i\) + ^iN{-iN), (6) 

and the sum of equations and (jS]) is identical to the parent equation ([2]). Random 
values {^°} and {^^} obey the identities = 2T_(5(ti - ta) and {^\ti)^\ti)) = 

2 AT5(ti — t2); is the total energy ([1]) where the arguments xi(t), X2(t), . . . , XAr(t) of the 
total process are replaced by coordinates of the equilibrium process x'j'(t), ^^(t), . . . ,x5^r(t). 
Expression in the square brackets in ([6]) is the difference of forces acting on the ith oscillator 
from the total process x(t) and equilibrium process x'^(t). It is worth mentioning that this 
force is the random value, and the process x^(t) (heat transport) is realized in the lattice 
with time- dependent random potentials. The problem of heat conductivity in the random 



time-independent potentials was analyzed in 



■421. 



Equation ([5]) describes the system embedded in the heat reservoir at temperature T_ 
of both lattice ends. And x°(t) is the stationary equilibrium process described by the 
canonical Gibbs distribution. Process x^(t) is responsible for the heat transport and the 
Wiener process {^^{t)} on the left lattice end defines temperature AT. Right lattice end 
has zero temperature (only the friction force acts on this oscillator). An expression for the 
local heat flux is 

J,^,+i = [F,^.+i(x) - F,^i+i(xO)] (7) 

and the equilibrium process x° does not transfer energy: {|Fj_i.j_,_i(x°) = 0, where (. . .) 

stands for the time average. It is essential that the heat flux ([71) is the small difference of 
large values from processes x(t) and x°(t). This is the reason why MD simulation gives large 
fluctuation when AT and T is not low. It can be shown that the time of computation 
increases oc (AT)~^ if the standard error is fixed. The comparison of two approaches (solving 
of standard SDEs (j2]) and ([5])- ([6])) is shown in Fig. [2] and results coincide with very good 
accuracy. 

Some results in this paper are obtained for the number of oscillators = 5. It may 
appear that this value is too small. For instance, "standard" simulations require up to 
~ 10'^ particles and ~ 10^ integration steps plus ensemble averaging But our results 
are aimed at establishing some new issues where number of particles is unessential. Lattices 
with larger number of oscillators were tested when necessary. 
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FIG. 2: Temperature dependence of heat conductivity for the lattice of = 5 oscillators. Filled 
circles: solution of standard SDKs ([2]); empty circles: SDKs (H])-®. Averaging over 100 MD 
trajectories 10^ time units (t.u.) each. T_ = 0.2, AT = O.OIT-. Triangle up at T = is the exact 
value in the harmonic approximation (/3 = 0). 

Relative displacements {x^ — x^_i) and {x} — x}_i) for processes x°(t) and x^(t) are shown 
in Fig. 3. These values characterize the energy fluxes in the lattice. Energy fluxes to the 
left and to the right are equal on average for process x'^(t) as it is the equilibrium process 
without energy transfer. But the energy flux is directed mainly to the right for process x^(t) 
(right panel). Low temperature is chosen for the better illustration. 

The dependence of heat conductivity on the oscillators number is shown in Fig.Hlat two 
value of temperature T_. Results coincide with very good accuracy. Inharmonicity becomes 
negligible at low temperature and heat conductivity at T_ = 0.1 (cirles in Fig. H]) coincides 
with the heat conductivity of the harmonic lattice (dashed line) with good accuracy. The 
analytical solution of the heat conductivity for the harmonic lattice is given in [43i]. There 
should be solved twice as large SDEs (I5])-(I6]) in suggested approach as that in standard 
scheme (E]). But this approach has some undoubted merits as discussed below. 

B. Heat conductivity at small temperature gradients 

One of the goals of the present paper is the computation of heat conductivity at very 
small temperature gradients. With this in mind an expression for the heat flux is analyzed 



7 





FIG. 3: Spatiotemporal evolutions of relative displacements {x^ — x^_i) and {xj — xj_i) for 
the equilibrium x^{t) (left panel) and nonequilibrium x^(t) (right panel) processes. N = 100, 
r_ = 0.05, AT = 0.0001. 
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FIG. 4: Coefficient of heat conductivity for the /3-FPU lattice for N = 7 — 150 oscillators. Squares: 
r_ = 1, circles: T_ = 0.1. Filled symbols - results obtained by the solution of standard SDEs 
([2]), empty symbols - SDEs ([S])-®. Averaging over 200 MD trajectories 310^ t.u. AT = 0.0ir__. 
Dashed line - harmonic approximation. Filled symbols are practically fully covered by empty 
symbols and are invisible. 

in more details. The expression for the local heat flux ([7]) can be rewritten as 

Ji^i+i = |'-^i^i+i(x) — Fj^j+i(x°)l x],-^ = 

(8) 

— Xi) — (Xj+i — Xi) + {Xj^_^i — Xj) + {x^_^_i — Xj) ] 
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where x}_^_i - velocity of {i + l)th oscillator in process x^. Taking in mind that Xi = + x}, 
the expression ([8]) can be transformed to 

J.^.+i = {-{xl, - xl) [1 + 3 {xl, - x^lf] + {xl, - xl f [1 + 3 (xl, - } xl,. (9) 

Processes x'^(t) and x^(t) have different ranges of specific energies. Noise terms {^^}, 



which provide temperature AT, are of the order ~ V AT (as (^^(ti)^^(t2)) ~ AT). And 
one can expect that process x^(t) has the same order x^(t) ~ V AT because equations (E]) 
become linear in the limit AT — when — )■ 0. Expression in curly brackets in is the 
polynomial of the third degree in the square root of temperature difference V AT. Taking into 



account that the velocity xj_^-^ is also of the order ~ vAT, expression iQ is the polynomial 



of the forth degree in y/ AT. But the coefficient of heat conductivity is determined by the 
relation J/ AT. Then terms of the third and forth orders can be neglected at AT — 0. Then 
is simplified to 

Ji^i+l = — (3^2+1 ~ ^i) 1 + 3 {x^_^_i — Xj) ii+i- (10) 

An expression for the potential energy, corresponding to process x^(t), can be derived 
analogously. This energy is the difference of potential energies f/(x) — [/'^(x°) and again, 
using coordinates x and x°, and preserving only terms quadratic in {x}_^_i — xj), one can get 
the potential energy for process x^ in the form 

= ^ E^^+i(^) (^'+1 - ^')'' 9^-,l{t) = 1 + 3 [xl,{t) - xM' , (11) 

i 

where gi{t) are time- dependent random coefficients of rigidities determined by process x°(t). 

It is illuminating to note that in the limit AT — t- the problem of heat conductivity in the 
/3-FPU lattice is reduced to the harmonic lattice with random coefficients. Corresponding 
SDEs have noise terms with friction forces on the left oscillator and zero temperature (only 
viscous forces) on the right oscillator: 



-9t{xl - Xi_i) + gi+i{x]_^i - x]) + 6ii{^^ - x\) - 6iNxlj (12) 



and QijQi^i are defined in (fTTj) . If the ID lattice with an arbitrary interaction potential is 
analyzed then the corresponding equations are the same with rigidities Qi = f/"(x° — x^_i) 
where U is the potential energy. SDK for the system with arbitrary neighbor radius of 
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interaction can be written in the general form as 

M 

i] = -Y, ^% ^] + -■*!)- ^^NX^N , (13) 

where A^^- - matrix of second derivatives of potential energy depending on x'', and M is the 
number of neighbors. 

C. Unusual dynamics of process x^(t) at high temperatures 

The process ^^(t) can be characterized by some time average correlators. The correlator 
([x\{t)]^^ was analyzed for the better understanding of the dynamics of process x^(t). One 
can expect that ([x\{t)]^^ ~ AT as discussed above. Two temperatures of the background 
process T_ were tested: Ti = 0.2 and T2 = 5 (hereafter subindex '-' is omitted, that is 
T_ = T). Results are shown in Fig. [51 

As expected, the correlator ( [xj(t)]^) linearly depends on AT: ( ~ AT at Ti = 
0.2. But the case is quite different at T2 = 5: ( reaches the stationary value ~ 0.064 
in the limit AT — > 0. It means that there exists some undamped stationary process x^(t) 
at high temperatures T of the background process x°(t) even in the limit AT 0. These 
results also imply an existence of a threshold temperature Tthr separating two regimes - 
damped at low temperatures and undamped at high temperatures. 

D. Threshold temperature 

Any process x^(t) damps out at low temperatures and flattens out to a stationary value 
at higher temperatures even in the limit AT — t- 0, and the temperature T of process x°(t) 
defines different damping rates. Bearing this in mind, it is convenient to excite some auxiliary 
process x^(t) over the background process x'^(t) and to analyze it. 

Coordinates and velocities of process x^(t) get random increments ^J^ilxjit = 0)]^ + 
^"^ilviit = 0)]^ = 0.5. The particular choice of initial conditions does not influence the 
final results. The total dynamics is the sum of two processes x(t) = x°(t) + x^(t). 

Stochastic differential equations for the process x^(t) are 

Xi = - 



dU dU" 



(14) 
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FIG. 5: Dependence of correlator ((a;})^) on the temperature difference AT. Filled circles: 
T = 0.2, empty circles: T = 5. Asymptotic value {{^i)'^) ^x^o ~ 0-064 at T = 5 (coefficient 
of linear regression 0.9993). Averaging over 100 MD trajectories 10'^ t.u. each. N = 5. The range 
of AT: 10-11 < AT/T < 2-10~i. 

and only viscous forces acts at the extreme left and right oscillators. U and are poten- 
tial energies with coordinates x(t) and x''(t), correspondingly. Stochastic dynamics f|T^ is 
imphcitly ruled out by the temperature T of process x°(t). 

To find the threshold temperature we initially consider the case of small temperature T 
when process Sl^{t) is damped out. Its damping is determined by the viscous friction of left 
{—Xi) and right (— x^) oscillators in ( |T^ . Gradually increasing the temperature its threshold 
value can be found when process x^(t) becomes undamped. The damping of mean squared 
displacement of the first oscillator [x}(t)]^ was calculated. Process x^(t) exponentially decays 
( [^i(^)]^) exp(— at) and a depends on T (see Fig. [6^1). The temperature dependence of 
coefficient a is shown in Fig. [6b and Tthr — 4.07 when a = 0. 

Next method to find the threshold temperature is moving "from up to down", going 
from higher to lower temperatures. At high temperatures there exists the stationary process 
arising from random forces = [dU/dxi — dU^ /dxf\ (see ([HD). And process x^ (t) decreases 
in a sense that all quadratic mean values tend to zero as temperatures approaches Tthr from 
above. When the temperature reaches its threshold value, process xi(t) is totally damped 
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FIG. 6: a) Exponential damping of process x^(t) at different temperatures: T = 3.5 (circles), 
T = 3.8 (squares), T = 4.0 (triangles up), T = 4.2 (triangles down). Solid lines - linear regressions. 
Averaging time ~5 000 — 10 000 t.u. 20 trajectories x'' were used to estimate the standard error, 
b) Damping coefficient (—a) as the function of temperature T of process x'^(t). Damping stops 
(a = 0) at Tthr ^ 4.07. N = 5. 
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FIG. 7: Stationary values ((a^x)^) at T > Tthr- Time of averaging 10^ t.u. The temperature 
dependence is approximated by the function ([2?i(t)]^) ~ exp[— 6/(T — Tthr)] (solid line). N = 5. 

(see Fig. [7]). The found threshold temperature is Tthr — 4.09. 

Strange behavior of process x^(t) is basically explained by time- dependent random forces 
$j = [dU /dxi — dU^ /dxf\ (see (E])) rather then random Langevin forces ~ V AT. And 
the plateau for the correlator ([x|(t)]^) ~ 0.064 at AT — )■ is determined exclusively by the 
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FIG. 8: Temperature dependence of the mean squared displacement ([a^i(i)]^)- Circles - MD 
simulation of SDEs (jl2p : dashed line - model of mean rigidities in the harmonic approximation. 
Averaging over 20 trajectories 2 10^ t.u. each. 

background process x°. 

This conjecture can be additionally supported. Let us consider the case of low tempera- 
tures T when process x° is "weak" . Then the rigidity coefficients Qi in ( ITT]) are close to unity. 
The lattice where actual rigidity coefficients gi (fTTl) are substituted by the mean values taken 
from the equilibrium Gibbs distribution: (gi) = go(T) and go(T) = 1 + 3 — is 
considered as an example. This harmonic model is exactly solvable and results are shown 
in Fig. E 

One can see that process damps out in the model with constant rigidity in 

contrast to the case when actual values (|TT1) are used. And the growth of process ([a^K^)]^); 
when temperature increases, is mainly governed by fluctuations rather then the increasing 
rigidities. 

Additional evidence of threshold phenomena is the one-dimensional analogue of the Math- 
ieu equation 

X = -[1 + gcos'^{t)]x. (15) 

Different types of solutions depend on the parameter g and initial conditions. There exists 
such critical value gcr that the solution is the superposition of periodic functions at 5^ < g^, 
and the solution diverges oc exp(±/xt) at g > gcr- 
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We consider an equation for the harmonic oscillator for one variable x with friction force 

X = —k(t)x — X, (16) 

where k{t) - stochastic rigidity. This equation is similar to equation f|T2|) for process x^(t) if 
k{t) = 1 + 3z'^{t) and z{t) is the stochastic process generated by the dynamics of harmonic 
oscillator with noise term and friction force at temperature T: 

z = -z + ^-z (17) 

and spectral property {^{ti) ^{t2)) = 2T6{ti — ^2)- The substitution a;exp(— 1/2) — )■ X 
excludes damping and (fT6l) becomes 

X = -[l + 3z\t)]X, (18) 

what is similar to the Mathieu's equation f[T^ . 

Eqs. (IT7|) - (IT5]) have rich family of solutions depending on initial conditions, temperature 
T and the sequence {^}. The solution is nearly harmonic function at very low temperatures. 
The superposition of harmonic functions is the solution at higher temperatures. At last 
there exists such threshold temperature Tthr that the solution diverges and is the product of 
harmonic functions by exp(/it). The solution is the product of some stochastic process by 
exp(z/t) at much higher temperatures, and u > fi. There are many interesting intermediate 
solutions, and this problem merits more attention. Considered examples show that an 
existence of threshold phenomena is not exceptional and can occur in different dynamical 
systems. 

The threshold temperature Tthr ~ 4.1 was found for the lattice length = 5. Larger 
lattice lengths were considered and the dependence of Tthr on the lattice length is shown 
in Fig. [9l Fitting gives dependence Tthr ~ 6 ■ 10^ A^"^. 

E. Time-resolved dynamics of process x^(t) 

Dynamics of process x^(t) at high temperatures T was analyzed above in terms of time- 
average correlators. And time- resolved behavior of process x^(t) at two temperatures T of 
the background process x° is shown in Fig. [TOl As above, A(t) = [2;J(t)]^ was calculated. 

One can see that A{t) behaves highly irregular. And numbers and heights of observed 
peaks increase with the growth of temperature until the process becomes chaotic at high T. 
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FIG. 9: Dependence of Tthr vs. lattice length N in log-log coordinates. Solid line is the fitting 
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FIG. 10: Dependence of A(t) = on time at different temperatures of process x'^(t). Left 

panel: T = 4.3; right panel: T = 7.0. N = 5, integration step h = 0.01. Time average (A(t))|^~g'^ 
are shown in horizontal solid lines. 

At temperatures T < Tthr process Si\(t) consists of individual rare peaks what can point to 
the possibihty that the energy can be transmitted by impulses. 

III. SOUND, SOLITONS AND BREATHERS IN THE /3-FPU LATTICE 



In accordance with the conjecture on the soliton contribution to the heat conductivity 
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, l35|, an attempt was made to shed some light on this problem. The spatiotemporal 
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FIG. 11: Correlator y^Q{t) y^o-^-mit + 20) as a function of the lattice coordinate n. /3-FPU lattice, 
N = 101, T = 2. Arrow points to n = 50 

correlator [nkit) yk+m(t + r)] was calculated, where yi{t) = Xi(t) — Xi^i(t) is the relative 
displacement of neighboring particles in time instant t. Solitons, being highly correlated 
displacements of particles, can leave a trace in correlation function. Time shift r = 20 was 
fixed and spatial correlation were calculated. Result is shown in Fig. [TTl Two peaks in the 
correlation function, shifted hj m ^ ±25, are visible. The velocity of their propagation is 
Vp = m/r ^ 1.25. 



A. Sound velocity in the /3-FPU lattice 

The temperature dependence of the sound velocity in the /3-FPU lattice was discovered 



about a decade ago 



44j |. The sound velocity was estimated fsnd ~ vT+~a, where a - 
parameter of renormalized frequencies depending on the temperature. Asymptotic value of 
ihe sound velocity in the high temperature limit fsnd ~ 1.22 T^/^ was derived recently in 



451]. If this formula apply to T = 2 then fsnd ~ 1-45 what differs from Vp ~ 1.25 found from 
correlation functions above. 

Below we derive more accurate expression for the sound velocity at low temperatures. 



It was shown 



461] that in nonlinear systems there exists a spectrum of frequencies which 



are proportional to the harmonic ones, according to a well defined law. Then the /3-FPU 
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potential can be rewritten as 



n{y) = ( 1 + ly') \y' (19) 



and an expression in brackets can be replaced by an effective rigidity fcgff 

u{y) = kj-y". (20) 

As a result the lattice becomes harmonic and it is necessary to find fcefr- It can be done in 
terms of mean field approximation (MFA). Mean value of potential energy is 

(«p(?/)) = fce4 (2/'> ' (21) 

where (?/^) is the mean of . According to the virial theorem, mean values of potential 
and kinetic energies are equal in the harmonic lattice, that is (up) = (uk). But the identity 
(wk) = T/2 holds for ID systems. Then 

fce4 = \ (22) 

The condition of self consistency of the MFA is 

1 + \ = ^eff (23) 

and it follows that (y^) = Tjk^^ and substitution of this expression into ( 123|) gives the 
self-consistent equation for kes- 

1 + T/{2k,s) = fceff (24) 

with the solution 

1 fl T 
= 2 + V 4 + 2 ■ 

Eq. defines the rigidity coefficient for the harmonic lattice with the renormalized spec- 
trum depending on temperature T. Thereby the temperature renormalized sound velocity 




Vsnd = ykes = \ - + \ - + (m = 1) (26) 



and fsnd = 1-27 for T = 2 what coincides with good accuracy with Vp ~ 1.25 found 
from correlation functions. The high temperature asymptotic of the sound velocity fsnd ^ 
0.84 T^^'^\rp^-^^- The temperature dependence of sound velocity for temperatures in the range 
< T < 10 is shown in Fig. [13] and very good agreement between analytical and "experi- 
mental" results is observed. 
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FIG. 12: The temperature dependence of sound velocity Vsnd- solid line - dependence ([26]) : empty 
circles - MD simulation. 

If the mean field approximation is applied to the lattices with cubic nonlinearity u{y) = 
± ^y^ then no dependence of the sound velocity on temperature is expected. Really, 
u{y) = 1(1 ± ly) y^ and = (l ± | {y)) = 1 as {y) = 1. 

The obtained results point to the fact that the nonlinearity can be also ignored in de- 
scribing the thermodynamics of the /3-FPU lattice, at least at T < 10, and the harmonic 
lattice with renormalized rigidity coefficients (!25|) is an adequate model. This conjecture 
was checked at different temperatures T = 1, 2, 5, 10 by the comparison of the total energies 
computed in MD simulation of the /3-FPU lattice and its harmonic model with renormalized 
rigidity coefficients. Very good coincidence of both energies testifies this hypothesis. 

There ex.ts an accurate expre=..on for the specific potential energy (n,ean potential 
energy of one oscillator) |47j] 

(f/p) 1 \K,/,{q) + Ks/M 



- 1 



q = 1/(8T) (27) 



N 8 L 2iri/2(g) 

derived from the thermodynamics of the /3-FPU lattice; K - modified Bessel functions. Spe- 
cific potential energy fl27|) and (f/p) /N computed in harmonic approximation with rigidity 
coefficients (!25|) coincidence with good accuracy. 
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B. Solitons and breathers in the /3-FPU lattice 



The discrete /3-FPU lattice can be reduced to the modified Korteweg - de Vries (mKdV) 
equation in the continuum long- wavelength approximation (see Appendix). The mKdV 
equation has solutions in the form of solitons and moving breathers |48j |. 

Solitons of compression and elongation has the form 



y{z,t) = ±^Bsech { B 



(28) 



where plus/minus signs stand for elongation/compression solitons; B - single parameter 
which simultaneously determines amplitude, width and velocity of soliton; y - local defor- 
mation of the lattice; z - soliton coordinate at time t. 
The two-parameteric breather solution is 

, ^, . o , cos $ - iB/a) sin $ tanhvE^ 
y{z, T) = -4/3 sech^^ , ^ f o) L ■ 2^ > 29 

1 + [p/a)^ sm <P sechW 

where 

^ = 2l3{z --ft) + ilj, $ = 2a{z - 5t) + (j), and 
7 = 4(3^2-/32), 5 = 4(^2-3/32) 

and a,/3 are free parameters (see Appendix for more details). It is necessary to make 
transformation from continuous variables to discrete variables y —> Xi ~ xi^i and z i m 
an attempt to use soliton and breather solutions on the discrete /3-FPU lattice. 

If soliton ( l28l) and breather ( l29i) solutions exist in the continuum limit , then the question 
arises: whether these moving localized excitations can be observed on discrete lattice? The 
answer is 'yes' and below the visualization method is suggested. 



Recall that the visualization method for standing discrete breathers is well known [491. |50| : 
boundaries with friction forces absorb thermal noise (phonons) and standing breathers can 
be easily seen. Other method is necessary to visualize the moving excitations. Let we 
have the thermolized lattice with oscillators at temperature T. If "cold" lattice (with 
zero velocities and displacements) is switched to the thermalized lattice then solitons and 
breathers should "run out" to the cold lattice and could be observed. Results are shown in 
Fig. [13] and solitons and breather are immediately seen. 

Solitons and breather (shown in inserts to Fig. [T^ move faster then the sound front. At 
the first glance it is inconsistent with the relation between velocities of sound and solitons. 



19 




FIG. 13: Solitons and breather running out of the initially thermalized lattice. A ~ breather, B - 
soliton of compression, C - pair of antisolitons (solitons of elongation). Arrow at n = 200 shows 
the border separating initially thermalized and "cold" parts of the lattice. Initial temperature of 
the left lattice part (1 < n < 200), T = 10. 

The sound velocity at T = 10 is v^nd ~ 1.67 what is larger then the maximal soliton velocity 
(fsoi)max ~ 1-3 (see Appendix). But the temperature of expanding thermal excitations 
gradually decreases and the sound velocity also decreases according to ( !26|) . And there 
comes a point when solitons, which have constant velocity, keep ahead the sound front. 

There are good grounds for believing that solitons and breathers do exist in the /3-FPU 
lattice. Very likely that the soliton contribution to the heat conductivity increases with the 
growth of temperature. Really, the temperature dependence of the soliton density obeys the 



relation n{T) ~ T^/^ for the Toda lattice at low temperatures 5l|- Conceivably the growth 
of the solitons density with temperature might be an inherent property of nonlinear systems. 
Our results on the soliton contribution to the heat conductivity are inconsistent with 



previous publication 



45| where the energy carriers are effective phonons rather than solitons 



The possibility of energy transfer by solitons was conjectured three decades ago (28|. Less 



20 



studied is the possibility of e nerg y transfer by breathers. One suggested mechanism is the 



Targeted Energy Transfer 



52 



53| when an efficient energy transfer can occur under a precise 



condition of nonhnear resonance between discrete breathers. Various aspects and possible 



applications of energy transfer by breathers are considered in 



50|. 



IV. CONCLUSIONS 



In conclusion we briefly summarize our results. A new method for the calculation of 
the heat conductivity is suggested. This is done by the decomposing of the total dynamics 
into two parts: equilibrium process x''(t) at equal temperatures T of both lattice ends, and 
nonequilibrium process x^(t) at temperature AT of one end and zero temperature of the 
other. This approach allows to extract and analyze the heat conductivity in an explicit form. 

The primary goal of the paper was to develop a method which would allow to decrease the 
computational time at small temperature gradients when fluctuations of the heat flux are 
usually too large. It was supposed that at small temperature gradients, when the harmonic 
approximation is valid and an expression for the heat flux has the form (fTOl) . an analytical 
averaging over random Langevin noise terms can be done. This approach is very efficient 
for the calculation of quadratic in (t) terms - the gain was thousand- fold. But formulaes 
are very complex for the linear terms and an efficient algorithm for their realization was not 
found yet. By this expedient the objective has not been met in full: the gain in computational 
time is obvious for small (A^ < 100) lattices, but decreases as the lattice length increases. 
Nevertheless we suppose that the further analysis of process x^(t) can be useful as it is 
responsible for the energy transfer. 

The threshold phenomena are familiar in microcanonical ensembles [l^. There exists two 
values of specific energy E. One separates dynamical regime and weak chaos, and higher E 
separates weak and strong chaos. It may be inferred that an existence of threshold phenom- 
ena is also a common occurrence in canonical ensembles. Really, a threshold temperature 
Tthr was found. The threshold temperature separates the different behavior of process x^(t): 
process x^(t) damps out at T < Tthr and reaches the stationary value at T > Tthr- It niay 
be conceived that the soliton and breather contributions to the heat conductivity increases 
with the growth of temperature if T > Tthr- Solitons and breathers can emerge from either 
thermal fluctuations or higher order phonon interactions. Additional experiments for nano- 
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sized systems at low temperatures when T < Tthv can reveal some new features omitted in 
the present work. 

The modified Korteweg - de Vries equation is derived in the continuum approximation 
for the /3-FPU lattice. mKdV has solutions in the form of compression/elongation solitons 
and breathers. The stability of these quasi-particles was checked in numerical experiments. 
Both types of excitations were directly visualized. 

On the other hand, it was found that the non-linear /3-FPU lattice can be reduced to the 
harmonic lattice with the temperature renormalized frequency spectrum. This reduction 
allows to reproduce adequately the heat conductivity and thermodynamics of the parent 
lattice. These two, mutually contradictory, properties of the /3-FPU lattice, - an inherent 
existence of solitons and breathers, and its "harmonic" behavior, seem to be very strange. 
Further analysis is necessary to solve this dilemma. 

It is likely that some fraction of total heat conductivity is conditioned by solitons and 
moving breathers. But their contribution to the energy transfer deserves further investiga- 
tion. 

The /3-FPU lattice is unique in the sense that in the continuum limit it has stable solutions 
in the form of solitons of compression and elongation. If the number of compression and 
elongation solitons is equal on average, then no macroscopic changes in the lattice length 
appear and an additional energy of deformation is negligible. It is an additional energetic 
factor favoring the solitons existence. 

The case is quite different for unsymmetrical potentials which in the lowest order of 
the Taylor expansions have the form u{y) = ||/^ ± |?/^. a-FPU, Toda, Morse, Lennard- 
Jones potentials are the examples. All these "cubic" potentials can be easily reduced to the 
ordinary KdV equation with the solution in the form of soliton of compression. And the large 
number of solitons is highly energetically unfavorable due to macroscopical compression of 
the lattice length. 



Appendix A: The modifief Korteweg — de Vries equation 
for the /3-FPU lattice. Solitons and breathers. 

An approximate solution for the soliton of compression in the /3-FPU lattice was ob- 



tained about two decades ago [5J]. Analogous soliton solution with the profile Qsnd 
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-'snd 



2 

snd 



{z Xi — Xi-i) was used recently [45|] and this 
solution is the function of a single parameter - sound velocity fgnd- 

Below we derive an equation for the continuum analogue of the /3-FPU lattice with more 
accurate soliton and breather solutions. The /3-FPU potential has general form 



(Al) 



where yi = x—Xi^i is the relative displacement of neighboring oscillators. The corresponding 
equations of motion are 



(A2) 



Starting from (lA2p . the continuum approximation can be derived supposing small devi- 
ations from equilibrium. The series expansion in terms of yi up to the forth order is: 



yi±i =yi± y'i + ^y- ± ^y-' + ^y-^ . 



(A3) 



Substitution of this expansion into ( lA2p gives the continuum equation: 

l/ = «(y" + ^/^)+3/3(y^yO'- 



(A4) 



Below we follow the well known reductive perturbation method (RPM) 55|, |56[ to get nec 



essary equation in partial derivatives. The receipt consists in introducing new variables 



y = e^/^ui + e^/^s + • • • ; 
^ = e'/\z-ct)- 
T = e^'H. 



(A5) 



Next the hierarchy of the expansions in terms of small parameter e should be used. If 
[5| is substituted in ( ]A4I) and terms of the order e^/^ and higher are neglected, then 



= (^{ui}^^ and c = y/a. c is the sound velocity in the harmonic approximation. 
Transformation (lASP means that the new coordinate system ^ moves with velocity c relative 
to the old coordinate system z. 

Equation with the accuracy of the order e^/^ is 



a 



12 



(A6) 
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and one can see that (1A6P reminds the well known modified KdV equation. Additional 
variables substitution ui = y^a/Q/3w and r = {24/ ^/a)T should be done to get the exact 
form of the mKdV equation: 

wt + Gw'^w^ + w^^^ = . (A7) 
^ and T are spatial and time variables. 



Equation (IA7p has two types of solutions j48|. The soliton solution is 



w{C, T) = ±B sech(S^ - B^T + 5) = ±B sech[5({ - B'^T) - ^ . (A8) 

Plus/minus signs are related to elongation/contraction solitons, correspondingly. Soliton 
(lASP is the one-parametric solution: free parameter B defines simultaneously amplitude 
(B), half- width (~ and velocity (-B^); defines the soliton coordinate at T = 0. 
Returning back to initial coordinates {z,t}, the soliton solution is 

y{z, t) = ±J—B sech {B [z - {l + ^B^ /2A) t] - zo] , (A9) 

where z - coordinate, t - time, zq - soliton coordinate at t = 0. The soliton velocity is 
Vso\ = v^(l + and is "supersonic" relative to the sound velocity in the harmonic 

approximation v^^^ = i/a. If B increases then soliton has larger amplitude, becomes more 
narrow and its velocity increases. The soliton solution for the /3-FPU lattice can be written 
if discrete variables in f lA9P are used: Hi = Xi — z ^ i, zq ^ iq. 

Parameter B can be arbitrary large in the continuum limit f lA9|) . But the lattice dis- 
creteness imposes limitations on the soliton width: solitons with the half-width less then 
< 2.0 become unstable. That is to say, soliton amplitude and velocity also have upper limit: 



A = y'a/QPB < 1, 1 < t>soi < 1.3 and the free parameter B for discrete /3-FPU lattice 



Two-parametric breather solution is 

cos$ - (/S/a) sin\&tanh\& 



w(^,T) = -4/3 sech ^ 
where 



1 + (/3/a)2sin^$ sech^ 



(AlO) 



^ = 2/3^ - 7T - ^, $ = 2a^ -5T-(t)] 

(All) 

7 = 8/3(2a2-/32), 5 = 8a(a2 - 3/3^) 
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FIG. 14: Movement and collision of two solitons and one breather. Snapshots are made at 
t = 0, 100, 200, 350 t.u. 

a, (3 are free parameters; ip, 4> - initial phases; the group and phase velocities are Vg^ = 7/2/3 
and fph = S/2a, correspondingly. Returning back to coordinate system z, the breather 
solution takes the form 

cos<l> - sin$ tanh\E' 



y{z,t) 




^ 13 sech^ 
3/3 1 + (/3/a)2sin2$ sech^ 



where 



m = 2(3 



z — Joi 1 



7 
48/3 



$ = 2a 



Z — JOL 1 



48a 



(A12) 



(A13) 



and 7, b are defined in (lAlip . 

The group breather velocity = \fa.{\ —7/48/3) = ^Ja [1 — (2a^ — /3^)/6]. Its amplitude 
~ v/8a/3/3. 

Soliton and breather stability was checked in numerical simulation. Initial conditions 
are chosen in the form of two solitons and one breather. Relative displacements are chosen 
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according to (1A9I) and ( lAlOl) . correspondingly. Velocities were obtained by differentiating 
by time. 

The most left soliton of elongation has amplitude Ai ~ 1.09 and velocity vi ~ 1.3 and 
initially located at ii = 20. Next soliton of compression has amplitude A2 — 0.63 and velocity 
V2 — 1.1 and initially located at ^2 = 50. Breather with parameters a = tt/6, f3 = 1/6 is 
initially centered at is = 100 and its velocity V3 = 0.91 (see Fig. 14a). The velocities 
condition vi > V2 > V3 says that solitons and breather should meet and interact. This triple 
interaction is shown in Fig. 14c. But as result, all three species preserve their individuality 
(Fig. 14d). These findings (large free paths, preserving amplitudes and shapes after collision) 
unambiguously demonstrate that soliton and breather are do really exist in the /3-FPU lattice 
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